Context

[…]

1. Information about the dataset

Link to dataset

“Features are computed from a digitized image of a fine needle aspirate (FNA) of a breast mass. They describe characteristics of the cell nuclei present in the image.”

fine needle aspirate (FNA) of a breast mass

\(~\) \(~\)

Ten real-valued features are computed for each cell nucleus:

The mean, standard error and “worst” or largest (mean of the three largest values) of these features were computed for each image, resulting in 30 features. All feature values are recorded with four significant digits.

The 3-dimensional space is that described in: [3].

This database is also available through the UW CS ftp server: ftp ftp.cs.wisc.edu cd math-prog/cpo-dataset/machine-learn/WDBC/

2. Question of interest

The aim is to predict whether the cancer is benign or malignant with only few features

3. Choice of data

There is a lot of data, we choose only the data concerning the mean. The following analysis has been done also for the worst-type of data and lead to similar conclusion. See in annex the inspection of these 2 subframes.

4. Plan of the analysis

• Inspect the data to understand it better
• Check if the sample proportion of benign and malignant cancer is representative of the whole population
• Select variables of interest by discarding according to the correlation level
• Built GLM models and select the best model with anova, discuss the goodness of this model
• Apply PCA to see to what extent dimension can be reduced and perform a GLM model with relevant dimension.
• Conclude by addressing the accuracy of these GLM models with a subset of the data

In annex, I built a GLM model with the library glmnet and I find the same selection of variable.

1. Required libraries

# Check if packages are installed and if not to install them or just load them
if(!require(pacman)) {
  install.packages(c("pacman", "remotes"))
}
if (!require(papaja)) {
  remotes::install_github("crsh/papaja")
}
if(!require(pacman)) {
  install.packages("equatiomatic")
}
if(!require(devtools)) {
  install.packages("devtools")
  install_github("kassambara/factoextra")
}
pacman::p_load(pander,       # Rmarkdown visualization
               GGally, 
               ggfortify,
               ggplot2,
               MASS,         # for stepAIC
               here,
               kableExtra,   # Rmarkdown table visualization
               papaja,
               glmnet,       # GLM implementation and analysis library
               equatiomatic, # Rmarkdown model equation visualization
               patchwork,    # arrange subplot
               devtools,     # tool for PCA
               factoextra,   # tool for PCA
               
)

## 
## Your package installed

## Warning in pacman::p_load(pander, GGally, ggfortify, ggplot2, MASS, here, : Failed to install/load:

2. Import data

path = here("LUCILE") # get relative path
setwd(path) # set working directory
df <-
  read.csv('data.csv', stringsAsFactors = 1)# load data from github repository

Let’s delete ID number and the last variable because there are not relevant.

df<-df[,-33]
df<-df[,-1]

We will work only with the mean-type data, so let’s create a new frame for the variables of this type.

4. Proportion of benign vs malignant cancer

kable(prop.table(table(df_mean$diagnosis)),col.names = c("Severity","Frequency"))   

The two types of cancer are not represented in the same proportion, this can lead to a bias.

However, this proportion representative of the reality:

“The benign to malignant ratio (B:M ratio) among breast biopsies (number of benign breast lesions divided by number of breast cancers) is widely believed to be around 4:1 or 5:1” [2]

2. Variables selection

• As expected, radius, perimeter and area are highly correlated (\(r \approx 1\))

• Surprisingly, concavity, compactness and concave.points have a strong correlation.( \(r \approx 0.8\) or \(r \approx 0.9\))

• Even if compactness is define as \(perimeter^2 / area - 1.0\)m the r between this variable and area or perimeter is not 1 because the correlation show only the linear dependency and their relation is not linear.

m_perimeter <-
  lm(
    data = df_mean,
    perimeter ~ radius + texture + area + smoothness + compactness + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_perimeter))

m_area <-
  lm(
    data = df_mean,
    area ~ radius + texture + perimeter + smoothness + compactness + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_area))

m_compactness <-
  lm(
    data = df_mean,
    compactness ~ radius + texture + perimeter + smoothness + area + concavity +
      concave.points + symmetry + fractal_dimension
  )
pander(summary(m_compactness))

df_mean_reduc <-df_mean[-c(4,5,7)]

1. Set a first GLM model

We set a GLM model with the remaining features. why, why binomial : […]

m <-
  glm(
    data = df_mean_reduc,
    diagnosis ~ radius + texture + smoothness + concavity + concave.points +
      symmetry + fractal_dimension,
    family = binomial
  )
summary(m)

## 
## Call:
## glm(formula = diagnosis ~ radius + texture + smoothness + concavity + 
##     concave.points + symmetry + fractal_dimension, family = binomial, 
##     data = df_mean_reduc)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.35180  -0.13938  -0.03229   0.02046   3.15368  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -28.38387    6.66946  -4.256 2.08e-05 ***
## radius               0.88701    0.21852   4.059 4.93e-05 ***
## texture              0.37262    0.06212   5.998 2.00e-09 ***
## smoothness          78.50170   32.64920   2.404   0.0162 *  
## concavity           15.52082    8.35462   1.858   0.0632 .  
## concave.points      46.67203   26.16265   1.784   0.0744 .  
## symmetry            16.85783   10.75613   1.567   0.1170    
## fractal_dimension -101.54448   61.26233  -1.658   0.0974 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 751.44  on 568  degrees of freedom
## Residual deviance: 153.35  on 561  degrees of freedom
## AIC: 169.35
## 
## Number of Fisher Scoring iterations: 8

3. Validation

The selected model is the following:

m_selected <-
  glm(
    data = df_mean_reduc,
    diagnosis ~ radius + texture + smoothness + concave.points,
    family = binomial
  )
pander(summary(m_selected))

(Dispersion parameter for binomial family taken to be 1 )

With the formula:

extract_eq(
  m_selected,
  use_coefs = TRUE,
  wrap = TRUE,
  terms_per_line = 2
)

\[ \begin{aligned} \log\left[ \frac { \widehat{P( \operatorname{diagnosis} = \operatorname{M} )} }{ 1 - \widehat{P( \operatorname{diagnosis} = \operatorname{M} )} } \right] &= -28.58 + 0.85(\operatorname{radius})\ + \\ &\quad 0.36(\operatorname{texture}) + 52.26(\operatorname{smoothness})\ + \\ &\quad 78.74(\operatorname{concave.points}) \end{aligned} \]

The ratio of the residual deviance by its degrees of freedom is \(160.3/564 = 0.284\) where the dispersion parameter is 1. [we want these two values close…]

diagnostic plot

autoplot(m_selected,1:4)

[…]

See in annex an other GLM modelisation with glmnet less conservative.

2. Observations in PC plans

autoplot(p,x = 1,y = 2, data = df_mean_reduc,colour = "diagnosis",
  loadings = TRUE,loadings.colour = "blue",loadings.label = TRUE)+
  theme_minimal()

By plotting the data in the plan of the first 2 principal components we see a clear separation between benign and malignant type of cancer. That mean that knowing the location in this plan of a new observation should allow to predict the severity of the cancer. Let’s try this approach by applying a GLM model to the first two PC.

Note: We tried to plot the data in all the possible pairs formed by the first four components. The clearest separation between the type of cancer occurs in the plan PC1-PC2

3. GLM model with PCA

df_mean_pca <- cbind(df_mean_reduc, p$x)
glm_pca <-glm(data= df_mean_pca, df_mean_pca$diagnosis ~PC1+PC2, family = binomial)

pander(summary(glm_pca))

autoplot(glm_pca,1:4)+ theme_minimal()

1.Features allowing to predict severity of cancer ?

[…] ## 2. test accuracy of the prediction with these feature

1. Correlation in the ‘standard deviation’ and ‘worst’ group

## 'data.frame':    569 obs. of  11 variables:
##  $ diagnosis         : Factor w/ 2 levels "B","M": 2 2 2 2 2 2 2 2 2 2 ...
##  $ radius            : num  1.095 0.543 0.746 0.496 0.757 ...
##  $ texture           : num  0.905 0.734 0.787 1.156 0.781 ...
##  $ perimeter         : num  8.59 3.4 4.58 3.44 5.44 ...
##  $ area              : num  153.4 74.1 94 27.2 94.4 ...
##  $ smoothness        : num  0.0064 0.00522 0.00615 0.00911 0.01149 ...
##  $ compactness       : num  0.049 0.0131 0.0401 0.0746 0.0246 ...
##  $ concavity         : num  0.0537 0.0186 0.0383 0.0566 0.0569 ...
##  $ concave.points    : num  0.0159 0.0134 0.0206 0.0187 0.0188 ...
##  $ symmetry          : num  0.03 0.0139 0.0225 0.0596 0.0176 ...
##  $ fractal_dimension.: num  0.00619 0.00353 0.00457 0.00921 0.00511 ...

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## 'data.frame':    569 obs. of  11 variables:
##  $ diagnosis         : Factor w/ 2 levels "B","M": 2 2 2 2 2 2 2 2 2 2 ...
##  $ radius            : num  25.4 25 23.6 14.9 22.5 ...
##  $ texture           : num  17.3 23.4 25.5 26.5 16.7 ...
##  $ perimeter         : num  184.6 158.8 152.5 98.9 152.2 ...
##  $ area              : num  2019 1956 1709 568 1575 ...
##  $ smoothness        : num  0.162 0.124 0.144 0.21 0.137 ...
##  $ compactness       : num  0.666 0.187 0.424 0.866 0.205 ...
##  $ concavity         : num  0.712 0.242 0.45 0.687 0.4 ...
##  $ concave.points    : num  0.265 0.186 0.243 0.258 0.163 ...
##  $ symmetry          : num  0.46 0.275 0.361 0.664 0.236 ...
##  $ fractal_dimension.: num  0.1189 0.089 0.0876 0.173 0.0768 ...

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We get the same results as for the mean group. ## 2. STEP AIC

m5<-stepAIC(m,trace = TRUE)

## Start:  AIC=169.35
## diagnosis ~ radius + texture + smoothness + concavity + concave.points + 
##     symmetry + fractal_dimension
## 
##                     Df Deviance    AIC
## <none>                   153.35 169.35
## - symmetry           1   155.78 169.78
## - fractal_dimension  1   156.25 170.25
## - concave.points     1   156.63 170.63
## - concavity          1   156.64 170.64
## - smoothness         1   159.85 173.85
## - radius             1   173.68 187.68
## - texture            1   203.11 217.11

m5$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## diagnosis ~ radius + texture + smoothness + concavity + concave.points + 
##     symmetry + fractal_dimension
## 
## Final Model:
## diagnosis ~ radius + texture + smoothness + concavity + concave.points + 
##     symmetry + fractal_dimension
## 
## 
##   Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                        561   153.3485 169.3485

3. glmnet

Let’s separate the dataset into train and test set

## 75% of the sample size
prop_train_test <- floor(0.75 * nrow(df_mean_reduc))

## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(df_mean_reduc)), size = prop_train_test)
#sample = sample.split(data$num, SplitRatio = .75)
#train = subset(data, sample == TRUE)
#test  = subset(data, sample == FALSE)
train <- df_mean_reduc[train_ind, ]
test <- df_mean_reduc[-train_ind, ]

x_train <- train[,-1]
y_train <- train$diagnosis
x_test <-  test[,-1]
y_test <- test$diagnosis

tol_length=length(levels(y_train))

# run glm : train and test with cross validation
cvfit<-cv.glmnet(as.matrix(x_train), y_train,family = "binomial", type.measure="class")
plot(cvfit)

cvfit$lambda.min

## [1] 0.01890567

assess<-assess.glmnet(cvfit,newx=as.matrix(x_test), newy=y_test, s='lambda.min')
confusion.glmnet(cvfit, newx =as.matrix(x_test), newy = y_test, s = 'lambda.min')

##          True
## Predicted  B  M Total
##     B     79 11    90
##     M      1 52    53
##     Total 80 63   143
## 
##  Percent Correct:  0.9161

as.numeric(1-assess$class)

## [1] 0.9160839

coef(cvfit, s="lambda.min")

## 8 x 1 sparse Matrix of class "dgCMatrix"
##                            s1
## (Intercept)       -11.9549638
## radius              0.3363175
## texture             0.1750409
## smoothness          .        
## concavity           .        
## concave.points     63.6306102
## symmetry            .        
## fractal_dimension   .

Let’s look with the first lambda:

cvfit$lambda.1se

## [1] 0.05773519

assess<-assess.glmnet(cvfit,newx=as.matrix(x_test), newy=y_test, s='lambda.1se')
confusion.glmnet(cvfit, newx =as.matrix(x_test), newy = y_test, s = 'lambda.1se')

##          True
## Predicted  B  M Total
##     B     80 13    93
##     M      0 50    50
##     Total 80 63   143
## 
##  Percent Correct:  0.9091

as.numeric(1-assess$class)

## [1] 0.9090909

coef(cvfit, s="lambda.1se")

## 8 x 1 sparse Matrix of class "dgCMatrix"
##                            s1
## (Intercept)       -6.91104517
## radius             0.20042878
## texture            0.07001007
## smoothness         .         
## concavity          .         
## concave.points    42.10066651
## symmetry           .         
## fractal_dimension  .

4. knn

k nearst ? pas interpretable

#install.packages('FNN')
#install.packages('crosstable')
library(crosstable)
library('FNN')
dat_pred <- FNN :: knn(train =x_train,test=  x_test,cl= y_train, k =10)
pander(table(x = y_test, y = dat_pred))

Version of R used

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] FNN_1.1.3            crosstable_0.3.2     factoextra_1.0.7.999
##  [4] patchwork_1.1.1      equatiomatic_0.3.0   glmnet_4.1-3        
##  [7] Matrix_1.3-4         kableExtra_1.3.4     here_1.0.1          
## [10] MASS_7.3-54          ggfortify_0.4.12     GGally_2.1.2        
## [13] ggplot2_3.3.5        pander_0.6.4         devtools_2.4.3      
## [16] usethis_2.1.5        papaja_0.1.0.9997    pacman_0.5.1        
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-2   ellipsis_0.3.2     rprojroot_2.0.2    flextable_0.6.10  
##  [5] base64enc_0.1-3    fs_1.5.2           rstudioapi_0.13    farver_2.1.0      
##  [9] remotes_2.4.2      ggrepel_0.9.1      fansi_0.5.0        xml2_1.3.3        
## [13] codetools_0.2-18   splines_4.1.1      cachem_1.0.6       knitr_1.36        
## [17] pkgload_1.2.4      jsonlite_1.7.2     shiny_1.7.1        compiler_4.1.1    
## [21] httr_1.4.2         backports_1.2.1    assertthat_0.2.1   fastmap_1.1.0     
## [25] cli_3.1.0          later_1.3.0        htmltools_0.5.2    prettyunits_1.1.1 
## [29] tools_4.1.1        gtable_0.3.0       glue_1.4.2         dplyr_1.0.7       
## [33] Rcpp_1.0.7         jquerylib_0.1.4    vctrs_0.3.8        svglite_2.0.0     
## [37] iterators_1.0.13   xfun_0.26          stringr_1.4.0      ps_1.6.0          
## [41] testthat_3.1.1     rvest_1.0.2        mime_0.12          lifecycle_1.0.1   
## [45] scales_1.1.1       promises_1.2.0.1   RColorBrewer_1.1-2 yaml_2.2.1        
## [49] memoise_2.0.1      gridExtra_2.3      gdtools_0.2.3      sass_0.4.0        
## [53] reshape_0.8.8      stringi_1.7.5      highr_0.9          desc_1.4.0        
## [57] foreach_1.5.1      checkmate_2.0.0    pkgbuild_1.3.1     zip_2.2.0         
## [61] shape_1.4.6        rlang_0.4.11       pkgconfig_2.0.3    systemfonts_1.0.3 
## [65] evaluate_0.14      lattice_0.20-44    purrr_0.3.4        labeling_0.4.2    
## [69] processx_3.5.2     tidyselect_1.1.1   plyr_1.8.6         magrittr_2.0.1    
## [73] bookdown_0.24      R6_2.5.1           generics_0.1.0     DBI_1.1.2         
## [77] pillar_1.6.3       withr_2.4.2        survival_3.2-11    tibble_3.1.5      
## [81] crayon_1.4.1       uuid_1.0-3         utf8_1.2.2         rmarkdown_2.11    
## [85] officer_0.4.1      grid_4.1.1         data.table_1.14.2  callr_3.7.0       
## [89] rmdformats_1.0.3   forcats_0.5.1      digest_0.6.28      webshot_0.5.2     
## [93] xtable_1.8-4       tidyr_1.1.4        httpuv_1.6.4       munsell_0.5.0     
## [97] viridisLite_0.4.0  bslib_0.3.1        sessioninfo_1.2.2